#!/usr/bin/perl -w
# gmaps2kml.cgi - Maps the geocode onto Google Earth's KML
# Hirotaka Ogawa (hirotaka.ogawa at gmail.com)

use strict;
use CGI;

# 直交座標経由の測地系変換
# Nowral   PXI07463@nifty.ne.jp
# 99/9/10
#
# 定数
my $pi  = 4 * atan2(1,1); # 円周率
my $rd  = $pi / 180;      # [ラジアン/度]

# データム諸元
# 変換元
# (Tokyo)
my $a = 6377397.155;
my $f = 1 / 299.152813;
my $e2 = 2*$f - $f*$f;

# 変換先
# (WGS 84)
my $a_ = 6378137;        # 赤道半径
my $f_ = 1 / 298.257223; # 扁平率
my $e2_ = 2*$f_ - $f_*$f_;  # 第1離心率

# 並行移動量 [m]
# e.g. $x_ = $x + $dx etc.
my $dx = -148;
my $dy = +507;
my $dz = +681;

my $q = CGI->new;
$q->charset('utf-8');

my $ll = $q->param('ll') || '';
my $mode = lc $q->param('mode') || '';

my ($lat, $lng);
if ($ll =~ /(\-?[0-9]+(\.[0-9]+)?),(\-?[0-9]+(\.[0-9]+)?)/) {
  ($lat, $lng) = ($1, $3);
  ($lat, $lng) = tky2wgs($lat, $lng) if $mode eq 'tky';
  print $q->header('application/vnd.google-earth.kml+xml');
  show_kml($lat, $lng,
	   ($q->param('t') || 'Pin It On The Earth!'),
	   ($q->param('r') || ''));
} else {
  print $q->header('text/html');
  show_html();
}

sub show_kml {
  my ($lat, $lng, $title, $url) = @_;
  my $desc = $url ? "<description><![CDATA[<a href=\"$url\">Original Map Location</a>]]></description>" : '';
  print <<KML;
<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://earth.google.com/kml/2.0">
<Placemark>
  $desc
  <name>$title</name>
  <styleUrl>root://styleMaps#default+nicon=0x304+hicon=0x314</styleUrl>
  <Point>
    <coordinates>$lng,$lat,0</coordinates>
  </Point>
</Placemark>
</kml>
KML
}

sub show_html {
  print <<HTML;
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    <title>gmaps2kml.cgi: Maps the geocode onto Google Earth</title>
    <style>h1 { font-size: 14px; } body, input { font-size: 12px; font-family: verdana,arial,sans-serif; }</style>
    <script src="http://maps.google.com/maps?file=api&v=2&key=ABQIAAAA2SGB7kYAwJ7vNIiGUvsYgBQhPSrtY3c-aq-G7Tm93zGQQxtL5xTnwm2_6ZCVnME0pZzKvB7i8npggw" type="text/javascript"></script>
  </head>
  <body>
    <h1>Maps the geocode onto Google Earth <small>(SEE ALSO <a href="http://as-is.net/blog/archives/001028.html">http://as-is.net/blog/archives/001028.html</a>)</small></h1>
    <div id="map" style="width:600px;height:400px;"></div>
    <form action="gmaps2kml.cgi">
    <ul>
    <li>Latitude,Longitude:<br />
	<input type="text" id="ll" name="ll" size="50" value="" /></li>
    <li>Datum:
	<input type="radio" name="mode" value="auto" checked="checked" />Auto-select
	<input type="radio" name="mode" value="wgs" />WGS
	<input type="radio" name="mode" value="tky" />TKY</li>
    <li>Title: (optional)<br />
	<input type="text" name="t" size="50" value="" /></li>
    <li>URL: (optional)<br />
	<input type="text" name="r" size="50" value="" /></li>
    </ul>
    <p><input type="submit" value="Generate KML" /></p>
    </form>
    <script type="text/javascript">
    //<![CDATA[
    var map = new GMap2(document.getElementById("map"));
    map.addControl(new GLargeMapControl());
    map.addControl(new GMapTypeControl());
    var ll = document.getElementById("ll");
    GEvent.addListener(map, "moveend", function() {
      ll.value = map.getCenter().toUrlValue();
    });
    map.setCenter(new GLatLng(35.658586, 139.745450), 16);
    //]]>
    </script>
    <p></p>
  </body>
</html>
HTML
}

sub tky2wgs {
  my ($b, $l) = @_;
  my $h = 0.0;
  my ($x, $y, $z) = &llh2xyz($b, $l, $h, $a, $e2);
  ($b, $l, $h) = &xyz2llh($x+$dx, $y+$dy, $z+$dz, $a_, $e2_);
  ($b, $l);
}

sub wgs2tky {
  my ($b, $l) = @_;
  my $h = 0.0;
  my ($x, $y, $z) = &llh2xyz($b, $l, $h, $a_, $e2_);
  ($b, $l, $h) = &xyz2llh($x-$dx, $y-$dy, $z-$dz, $a, $e2);
  ($b, $l);
}

sub llh2xyz { # 楕円体座標 -> 直交座標
  my($b, $l, $h, $a, $e2) = @_;
  my($sb, $cb, $rn, $x, $y, $z);

  $b *= $rd;
  $l *= $rd;
  $sb = sin($b);
  $cb = cos($b);
  $rn = $a / sqrt(1-$e2*$sb*$sb);

  $x = ($rn+$h) * $cb * cos($l);
  $y = ($rn+$h) * $cb * sin($l);
  $z = ($rn*(1-$e2)+$h) * $sb;

  ($x, $y, $z);
}

sub xyz2llh { # 直交座標 -> 楕円体座標
  my($x, $y, $z, $a, $e2) = @_;
  my($bda, $p, $t, $st, $ct, $b, $l, $sb, $rn, $h);
  $bda = sqrt(1-$e2); # b/a

  $p = sqrt($x*$x+$y*$y);
  $t = atan2($z, $p*$bda);
  $st = sin($t);
  $ct = cos($t);
  $b = atan2($z+$e2*$a/$bda*$st*$st*$st, $p-$e2*$a*$ct*$ct*$ct);
  $l = atan2($y, $x);

  $sb = sin($b);
  $rn = $a / sqrt(1-$e2*$sb*$sb);
  $h = $p/cos($b) - $rn;

  ($b/$rd, $l/$rd, $h);
}

sub deg2dms {
  my($d) = @_;
  my($m, $s, $sf);
  $sf = int($d*360000 + 0.5);
  $s = $sf / 100 % 60;
  $m = $sf / 6000 % 60;
  $d = int($sf/360000);
  $sf %= 100;
  sprintf("%d\.%d\.%d\.%d", $d, $m, $s, $sf);
}
